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Abstract We consider the Clairaut theory of the equilibrium ellipsoidal figures for differentiated non- 
homogeneous bodies in non-synclrronous rotation (Tisserand, Mecanique Celeste, t.II, Chap. 13 and 
14) adding to it a tidal deformation due to the presence of an external gravitational force. We assume 
that the body is a fluid formed by n homogeneous layers of ellipsoidal shape and we calculate the 
external polar flattenings e*,, /x& and the mean radius Rk of each layer, or, equivalently, their semiaxes 
dk,bk and Ck■ To first order in the flattenings, the general solution can be written as ek = Rk^h and 
/ifc = RklJ-hi where Rk is a characteristic coefficient for each layer which only depends on the internal 
structure of the body and eh , Rn are the flattenings of the equivalent homogeneous problem. For the 
continuous case, we study the Clairaut differential equation for the flattening profile, using the Radau 
transformation to find the boundary conditions when the tidal potential is added. Finally, the theory 
is applied to several examples: i) a body composed of two homogeneous layers; ii) bodies with simple 
polynomial density distribution laws and iii) bodies following a polytropic pressure-density law. 

Keywords Polar flattenings • Tidal potential • Rotation • Differentiated bodies • Clairaut equation • 
Ellipsoidal figure of equilibrium ■ Exoplanets • Polytropes 


1 Introduction 

Several theories of tidal evolution, since the theory developed by Darwin in the XIX century (Darwin, 
1880), are based on the figure of equilibrium of an inviscid tidally deformed body (see e.g. Ferraz-Mello 
et al., 2008; Ferraz-Mello, 2013). The addition of the viscosity to the model is done at a later stage, 
but the way it is introduced is not unique and can vary when different tidal theories are considered. 
Frequently, the adopted figure is a Jeans prolate spheroid or, if the rotation is important, a Roche 
triaxial ellipsoid (Chandrasekhar 1969). It is worth recalling that ellipsoidal figures, are excellent first 
approximations, but not exact figures of equilibrium (Poincare, 1902; Lyapunov, 1925, 1927). Besides, 
Maclaurin, Jacobi, Roche and Jeans ellipsoids are valid only for homogeneous bodies. Real celestial 
objects, however, are quite far from being homogeneous. This causes significant deviations which need 
to be taken into account in the astronomical applications. 

The non-homogeneous problem, when we only consider the deformation by rotation, has been 
extensively studied. The problem of one body formed by n rotating homogeneous spheroidal layers 
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as well as its extension to the continuous case was studied by Clairaut (1743) (revisited by Tisserand 
(1891) and Wavre (1932)). Their works were based on the hypotheses of small deformations (linear 
theory for the polar flattenings) and constant angular velocity inside the body. The general case 
of homogeneous layers rotating with different angular velocities (non-linear theory) was studied by 
Montalvo et al. (1983) and Esteban and Vazquez (2001) (see Borisov et al. 2009 for a more detailed 
review), and was generalized to the continuous inviscid case by Bizyaev et al. (2014). 

The case of uniformly rotating layers was studied by several authors. Kong et al. (2010) discussed 
the particular case of a body formed by two homogeneous layers with same angular velocity. Hubbard 
(2013), with a recursive numerical form of the potential of a N-layers rotating planet, in hydrostatic 
equilibrium, showed a solution for the spheroidal shapes of the interfaces of the layers. 

When the tidal forces acting on the body are taken into account along with the rotation, the 
literature is much less extensive. Usually the spin-orbit synchronism is assumed, so that the rotating 
body solution can be used (e.g. Van Hoolst et al. (2008)). Tricarico (2014), assuming synchronism, 
found a recursive analytical solution for the shape of a body formed by an arbitrary number of layers. 
For this, he developed the potentials of homogeneous ellipsoids in terms of the polar and equatorial 
shape eccentricities. However, the results do not include tidally deformed bodies whose rotation is 
non-synchronous, as, for instance, the Earth, solar type stars hosting close-in planets and hot Jupiters 
in highly eccentric orbits. 

In this work, we generalize the linear Clairaut theory, adding a tidal potential due to the presence of 
an external body, without the synchronism hypothesis. The layout of the paper is as follows: In Section 
2 we present the 2n classical equations of equilibrium. The resolution of the system of equilibrium 
equations is shown in Section 3. In Section 4, we study the Clairaut’s equation for the continuous 
problem an its solution. In Section 5 we calculate the potential at a point in the space due to the 
deformed body and we calculate a generalized Love number for the differentiated non-homogeneous 
bodies. In Section 6, we apply the theory to a body composed of two homogeneous layers, while 
bodies with continuous density laws are studied in Section 7. Finally, in the Section 8 we present the 
conclusions. 


2 Equilibrium equation of a fluid in rotation 

We consider a rotating inviscid fluid of mass uit and a mass point M orbiting at a distance r from 
the center of the primary in a plane perpendicular to the rotation axis. We assume that the fluid is 
composed of n homogeneous layers of density pj~ (k = 1, • • •, n), that each layer has an ellipsoidal shape 
with external semiaxes a*,, b & and c& along the coordinate axes, and angular velocity $7. We define 
the mean radius of each layer as Rk = yfaubkCk- We choose a reference system such that r = rx and 
ft = f2z where x, z are unit vectors along the axes x, z (Figure [1]). 



Fig. 1 Body of mass mr formed of n homogeneous layers of density pk and mean radius Rk rotating with 
angular velocity fl = Qz and a point mass M orbiting at a distance r from its center in a plane perpendicular 
to the rotation axis. 









3 


Now, if we consider one point on the surface of the £-th layer, with position vector = iex+rjty+Oz 
and velocity is = ft x x^, we can use the same equation used in the study of equilibrium ellipsoids (see 
Tisserand, 1891, Chap. 8 and 13; Jeans, 1929, Sec. 215-216; Jardetzky, 1958; Chandrasekhar, 1969), 
which expresses the fact that the total force acting on a point of its surface must be perpendicular to 
the surface 

fJ t$t oc V^Vg + ft x (ft x x^) (1) 

where 

Mb,vt,Ct) = Q + § + 4-i = o ( 2 ) 

a t °t L t 

is the equation of the surface of the ellipsoid, Vrj is the potential of the gravitational forces at x^ and 
the last term corresponds to the centripetal acceleration. The use of the above equilibrium equation in 
a case where the tidal force field is changing because of the external body needs a justification. Eqn. 
(1) means that no change in the shape of the body occurs because of internal forces; the shape will 
change, but only because of the relative change of the position of the external body. 

Hence, we obtain the equilibrium equations 
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for each layer. For this, we will use the 2 n equilibrium equations ©. 

The gravitational potential can be written as the sum of the potential due to the mass M and the 
sum of the potentials of each layer. It can also be written as the sum of the potentials of n superposed 
homogeneous ellipsoids, with semiaxes a*,, b k and c k and densities 

&k — Pk (6) 


with p n +i = 0. The mass of each partial ellipsoid is 

47T q 

m k = (7) 

If we call 14 the potential of each ellipsoid, the total potential is 

n 

V G = V tid +^2 v k . (8) 

k= 1 

As the equilibrium equations ([3|) are linear in V, we can write 
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3 Flattenings of the layers 


The next step is to calculate the contribution of each potential to the equilibrium equations If 
we consider the contributions to the potentials due to the inner and outer layers on the £-th layer, we 
obtain the equations 
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(see Appendix A in the Online Supplement), which can be solved with respect to the £-th layer 
flattenings, giving 
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where cm and ej are the flattenings of the equivalent MacLaurin and Jeans homogeneous spheroids: 


cm = 


5 R 3 n I2 2 
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ej = 


15 MR? n 
Aititt 3 


(13) 


where mx is the mass of the body. 

These flattenings are obtained as solutions of the rotational and tidal problem, respectively, when 
the deformed body is a homogeneous spheroid with the same mass mj- and mean radius R n as the 
considered body (see online Appendix B). The coefficients ae k , /3e k and 7 ^ are 
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If we divide the equations for ee by ej + £m and divide the equations for m by cm, we obtain the 
same equation for the two polar flattenings: 


where 


Hi = 


Rt_ 

Rn 


def 67 


3 e-i 


fe =1 


_ = 

+ cm cm 


+ ^2 PekH k , 

(15) 

k=e+i 


(.j = 1 

(16) 


It is worth emphasizing that these equations naturally associate the flattening of the homogeneous 
MacLaurin spheroid with the polar flattenings / i k calculated using the minor semi-axis of the tidally 
deformed equator. This is so because the tide also acts shortening the polar axis. While the flattenings e k 
increase because of the tide, the flattenings [i k remain the same as in absence of tide. Therefore, the tide 
increases the mean polar flattening of the layers. The 3 axes of the layer are a,j = Rj [1 +Hj (e«+2e j)/3]; 
bj = Rj [1 + Hj(e M — £j)/ 3]; Cj = i?j[l + Hj(—2e M — ej)/3]. 
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It is important to note that if the orbital motion is synchronous with the rotation, in which case 
ej = 3cm, the system m is completely equivalent to that found by Tricarico (2014), where the square 
of the polar and equatorial “ eccentricities ” used there are related to the polar flattenings through 
e pi « 2 ei and e 2 qi r* 2a - 2 /Xj. 

The calculations done are valid only for small flattenings, i.e. they assume that the perturbation 
due to the tide and the rotation are small enough so as not to deform too much the body (in the 
second order, the figure ceases to be an ellipsoid). 


4 Extension to the continuous case 

In order to extend to the continuous case (following Tisserand, 1891, Chap. l43h we assume that the 
number of layers tends to infinity so that the increments ARk = Rk — Rk-i are infinitesimal quantities. 
When ARk —> 0, the equation (fl5l) becomes 


r 2 O f rz=x rZ= 1 

f( x )n(x) = ^x 5 + p(z)d{z 5 n(z)) + x 5 / p(z)dR(z), (17) 

^ & Jz=0 Jz=x 


where x = R/R n is the normalized mean radius (:r( 0 ) = 0 in the center and x{R n ) = 1 on the surface), 
p(x) = p(R)/po is the normalized density (po is the density in the center, therefore p(0) = 1) and the 
function f(x ) is 


f{x) = 3 p(z)z 2 dz, 


(18) 


with /( 0 ) = 0 and /( 1 ) = /„. 

Deriving El with respect to we have 


r z= 1 
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(19) 


and deriving once more we obtain the differential equation for the flattening profile 
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( 20 ) 


It is a homogeneous linear differential equation of second order with non constant coefficients and 
it turns out to be the same for both flattenings. It is the same expression found by Clairaut (Jeffreys, 
1953). 

The Eq. (1171) allows us to calculate easily the limits R n that the proportionality parameter R can 
take at the surface. In the homogeneous case p(x) = 1, the integrals can be calculated trivially. At 
the surface x = 1, we obtain R n = 1. In the non-homogeneous case, if the density is a non-increasing 
function (|| ^ 0 ), we have, at the surface 
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( 21 ) 


Then, under the assumption of equilibrium, a non-lromogeneous body will have flattenings on the 
surface with values between 0.4 and 1 times the values they would have if the body was homogeneous. 


1 see Appendix o in the Online Supplement for more details. 
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4.1 Boundary conditions. Radau transformation 


The differential equation (1201) requires two boundary conditions to be solved. However, before attempt¬ 
ing to find these boundary conditions, we will show two relationships that will prove useful later. The 
first relationship is obtained from equation (1191) . where at x = 1 we have 

H' n = 2{l-'H n ). (22) 

The second relationship is obtained from the differential equation m by just multiplying it by 
f(x)/x 2 and evaluating the resulting equation in the neighborhood of x = 0 (Note that f(x) ~ 
x 3 + 3(dp / dx)$x A / A and p(x) ~ 1 + (dp/dx) qx). It is 


n’ = 


dpo Bq 

dx 4 


(23) 


where is the derivative of the density at x = 0. 

In practical applications (see Section 6), it is convenient to introduce the Radau transformation 


p{x) = 


x'H’(x) 


n(x) ’ 

and rewritten Clairaut’s equation as the Ricatti differential equation 


V 


V 


q(x) 


r] + q{x) = 0, 


where 


q(x) = 


def 6 / p{x)x 


f{x) 


- 1 


In the new variables the boundary condition is 

p(x = 0) = 0. 


(24) 


(25) 

(26) 


(27) 


The variable p is sometimes referred to as Radau’s parameter (Bullen, 1975). Defining p(x = 1) = p n 
and using the relationship (1!E1) and the transformation EH) . the boundary conditions of EU1) are 


n 


n 


2 

2 + Pn 
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2pn 
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(28) 


As a result of this relationship, if considering that 0.4 < < 1, we recover the classical result 

0 < p n < 3 (Tisserand, 1891). 

Finally, it should be noted that once p(x) is found, we may find the profile flattening from equation 
EH) . whose solution is 

H{x) = U n eK ^ z)/z dz . (29) 


5 Potential of the tidally deformed body 


The contribution of the fc-th ellipsoid to the potential at an external point x = xx + yy + zz is given 

(3 COS 2 IP, - 1) - (3cos 2 Ip 2 - 1), (30) 

where r* = |x| and and ^ are the angles between the direction of the point where the potential is 
taken and the coordinate axes x and y respectively. @ The total potential is the sum of the potentials 
of all ellipsoids: 
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2 For the details of the calculation of SVf, see Eq. (I A. 131) in Appendix A.l (in the Online Supplement). 
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where eh, Ph are the flattenings of the equivalent homogeneous ellipsoid and the constant kf is often 
called fluid Love number (Munk and MacDonald, 1960; Correia and Rodriguez, 2013). For a non- 
homogeneous body, we find 


, def 3 Xyfc=i m kR%Rk 
7 = 2 > 

or using the continuous model, 

kf = VT f P{z)d{z 5 n(z)). 

z Jn Jz —0 

Using the integral form of Clairaut’s equation to evaluate the integral, we have 


(32) 


(33) 


k f = 



(34) 


which shows the link of the fluid Love number with the coefficient R n . This relationship is based on 
the fact that both constants depend solely on the internal structure, characterizing the inhomogeneity 
of the body. In the homogeneous case 'H n = 1 thus recovering the classical result kf = 1.5. 


6 Two-layer Core-Shell model 


In this section we consider the simple case of a body formed of two homogeneous layers: a core with 
density p\ and mean outer radius R±, and a shell with density p 2 = \pi (with A < 1) and mean outer 
radius R 2 (Fig. [2]). The densities of the superposed ellipsoids are 
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Fig. 2 Density profile of a body formed by two homogeneous layers. £ is the mean outer radius of the core 
relative to the mean outer radius of the shell f? 2 - A is the shell density relative to the core density pi. 




Fig. [3ja 


Fig.[3jb 


Fig. 3 Possible values of Hi (core) and H 2 (shell) as functions of the core size £ and of the relative density of 
the shell A. 


Hence 

H)(a + (1 — A)£ 3 ) 2 

(2 + 3A) (2A + 5(1 - A)£ 3 ) - 9A(1 - A)£ 5 

2(A + (1 — X)e) (2 + 3A + 3(1 - A)£ 5 ) 

—A- ^(40) 

(2 + 3A) (2A + 5(1 - A)C 3 ) - 9A(1 - A)£ 5 

Fig. © shows the results obtained for the constants 'Hi and H.- 2 - We see that: 

— If A = 1 or £ = 1 the constants are Hi = H 2 solutions for a homogeneous body. 

— When the core is denser than the mantle, H 2 ^ Hi and the flattenings of the nucleus are smaller 
than the flattenings of the surface (where e\ = H\{ej + cm) ^ £2 = ^(ej + cm) and fi 1 = HiCm ^ 
H 2 = H 2 eM)- 

— Since H .2 ^ 1, the maximum surface flattening is given by the homogeneous solution. In presence 
of a core, the surface is always less flattened than it is in the homogeneous case. 

— While Hi may take all possible values between 0 and 1, H 2 is always larger than the critical limit 
0.4, corresponding to the degenerate limit case in which the whole mass would tend to concentrate 
in the center and would be surrounded by a zero-density shell (case of Huygens-Roche). Therefore 
the flattenings of the outer surface can never be less than 40% of the homogeneous reference values. 
This is the same result given by Eq. CD for the continuous case. 


Hi = 

H 2 = 
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Fig. 4 Possible values of kf as functions of the core size £ and of relative density of the shell A. 


6.1 Fluid Love number 

Using equation HMD , together with the expression for 77 2 (Eqn. HOI) , the expression of the fluid Love 
number kf is 

5(A + (1 - A)£ 3 ) (2 + 3A + 3(1 - A)£ 5 ) 

k f = ^-- -L- - l. (41) 

(2 + 3AJ f 2A + 5(1 - A)£ 3 J - 9A(1 - A)£ 5 

Figure HI shows the possible value of kf as a function of the core size £ and of the relative density 
of the shell A. If we obtain kf, for example by determining 77 n by direct observation of the surface 
flattenings, then equation (HD defines a continuous curve of possible values for the size of the nucleus 
£ and the relative density of the shell A under the hypothesis of two homogeneous layers. Moreover, as 
can be seen in this figure, it allows us to predict a maximum value for these physical parameters. 


7 Application to different density distribution laws 

In this section, we present some applications of the theory developed in this paper to bodies with 
continuous density distributions. For this we use two examples of density distributions: polynomial 
and polytropic density laws. 

In both cases the Clairaut’s equation is solved numerically after introduction of the variable defined 
by the Eq. m- The flattening profile 7~L(x) and the Love number are then obtained through the inverse 
transformation. 


7.1 Polynomial density functions 

We consider initially a simple polynomial density law: 


p{x) = l-x a , (42) 

where a > 0. Figure ([5] a) shows the density functions for a = 0.1,1, 2,10 and 100 as functions of the 
normalized mean radius x. 

The resulting flattening profiles T~L(x) are shown in Figure d5]b) . In all cases, the flattening profile 
'H(x) is an increasing monotonic function and for all x, the values of 77(x) increase when the power a 
increases. 

Note that, as discussed in Section 3, the value of 77„ is always greater than the limit value 0.4 
and less than 1. Particularly 77 n tends to 0.570 when a tends to 0 and 77„ tends to 1 when a tends 
to 00 (homogeneous case). The fluid Love number increases from 0.424 (when a tends to 0) to 1.5 
(when a tends to 00 ). These results can be seen in Figure where we also show the values of 
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Fig. 5 (a): Density profiles for polynomial density distributions with different values of a. (b): Flattening 
profile H(x) for the same density laws, a = 0.1 (black), a = 1 (red), a = 2 (green), a = 10 (blue) and a = 100 
(magenta). 


the flattening factor at the surface and the dimensionless moment of inertia C/ttitR This last 
parameter increases from 0.24 (when a tends to 0) to 0.4 (when a tends to oojfl- 



Fig. 6 Values of T-L„ (black), kf (red) and C/mr Ri, (blue) for different values of the exponent of the polynomial 
density law. 


7.2 Polytropics pressure-density laws 


We may consider a self-gravitating body in hydrostatic equilibrium with a more general polytropic 
pressure-density law: 

P = Kp 1+ $, (43) 

where P is the pressure, n is the polytropic index and I\ is constant. The differential equation for the 
density is then given by the Lane-Emden equation (Chandrasekhar, 1939) 


1 d ( 

dt) 


+ 0 n = O, 


3 


An elementary calculation allows one to find the relationship 


c 

mj' R ‘ 2 


2 In P^d-z 2 v 3+a 

3 p pz^dz 5 A 5+q- 


(44) 
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1 _J 

where p = 9 n and R = at; with a 2 = (n + 1 )Kpft /47 tG. The standard boundary conditions are 
0(0) = 1 and 0'(O) = 0. If 0 < n < 5 the solution 0(£) decreases monotonically and has a zero at a 
finite value £ = £i. This radius corresponds to the surface of the body where P = p = 0. 

It is worth mentioning that several real cases exist that correspond to polytropes. For example, when 
convection is established in the interior of a star the resulting configuration is a polytrope; when the 
gas is degenerate, the corresponding equations of state have the same form as the polytropic equation 
of state, etc. (see Collins, 1989). We also mention recent results by Leconte et al. (2011) showing that 
the density profile of hot Jupiters is well approximated by a polytrope. 




Fig. 0a 


Fig. 0b 


Fig. 7 (a): Density profiles for different values of the polytropic index, (b): Flattening profile R(x) for these 
density laws, n = 0.5 (black), n = 1 (red), n = 1.5 (green), n = 3 (blue) and n = 4.5 (magenta). 


Figure 0a) shows the density functions for n = 0.5,1.0,1.5, 3.0 and 4.5 as functions of the nor¬ 
malized mean radius x = R/a£ i obtained from the integration of the Lane-Emden equation. 

The resulting flattening profiles R(x ) are shown in Figure (0b). In all cases, the flattening profile 
R(x) is an increasing monotonic function and for all x, the values of R(x) decrease when the polytropic 
index n increases. 

As mentioned previously, the value of R n is always greater than the limit value 0.4. Particularly 
R n —> 0.4 when n —> 5. The fluid Love number decreases from 1.5 for n = 0 (constant density) to 0 
when n tends to the limit n = 5. These results can be seen in Figure 0), where we also show the values 
of the flattening factor and the dimensionless moment of inertia C/uitR n for values of n below 
the limit n = 5. The adimensional moment of inertia decreases from 0.4 (when n = 0) and tends to 0 
when n —>■ 5. 



Fig. 8 Values of Rn (black), fc/ (red) and C/mrR n (blue) for different polytropic indices n < 5. 
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8 Conclusions 

In this paper, we have extended the classical results on non-homogeneous rotating figures of equilibrium 
to the case in which the body is also under the action of a tidal potential due to the presence of an 
external body, without the restrictive hypothesis of spin-orbit synchronization. The only assumptions 
in this paper are a body formed by n homogeneous ellipsoidal layers in equilibrium and small enough 
tidal and rotational deformations with symmetry axes perpendicular to each other (remember that, in 
the second order, the figure ceases to be an ellipsoid). We have calculated the 2 n equilibrium equations 
for small flattenings and we have found that the two polar flattenings tk and pk are linearly related, 
both being proportional to the homogeneous reference values with a factor of proportionality Rk which 
is the same in both cases. The deformations propagate towards the interior of the body in the same 
way depending, in the first approximation, only on the density profile; it does not depend on the origin 
of the two considered deformations. Then the problem of finding the 2 n flattenings, correspond to 
finding the n coefficients T~ik with n equilibrium equations. An important consequence of this approach 
is that the flattening profile Hk is the same no matter if the rotation of the body is synchronous or 
non-synchronous and the results for Rk are the same found by Tricarico (2014). 

We have also studied the continuous case as the limit for a very large number of layers of infinitesimal 
thickness, which leads to the Clairaut’s differential equation for the function R{x) (i.e. the same 
equation for both flattenings). This result was expected because the coefficients of the Clairaut equation 
only depend on the internal distribution of matter p(x). Therefore, the differential equation that 
generates the functional form of the profile flattening R{x) does not change when we change the 
nature of the deformation, provided that it is small. For densities decreasing monotonically with the 
radius, we have found that, at the surface, R n takes values larger than 0.4 (see Eq. (ELD) and takes 
the limit value 1 in the homogeneous case. This means that the surface flattenings of a differentiated 
body are always smaller than the flattening of the corresponding homogeneous ellipsoids, but always 
larger than 40% of it. 

The results were applied to several examples. In the case of a body composed of two homogeneous 
layers the following results were obtained: 

— In a realistic case where the core is denser than the shell, the flattening of the nucleus is smaller than 
the flattening of the surface. This is a result classically known to Tisserand (1891) and discussed 
in recent papers by Zharkhov and Trubitsyn (1978), Hubbard (2013), Tricarico (2014). 

— In the presence of a core, the surface is always less flattened than the homogeneous reference 
flattening, but larger than 40% of the latter value. 

— The fluid Love number kf < 1.5, define a continuous curve of possible values for the size of the 
nucleus £ and the relative density of the shell A and predicts their maximum value. 

Finally, we studied bodies with different continuous density laws, first for some simple polynomial 
functions and then for polytropic profiles. The following results were obtained: 

— In all cases, the function R{x) is an increasing monotonic function. 

— For all x , the values of R{x) increase from 0.530 to 1 when the power a increases from 0 to oo, in 
contrast with the polytropic densities, in which the values of R[x) decrease from 1 to 0.4 when the 
polytropic index n increases from 0 to the limit case n = 5. 

— The fluid Love number kf varies between 0.326 and 1.5 in the same range of the power a for 
polynomial densities. For the polytropic laws, the fluid Love number kf varies between from 1.5 to 
0 when the polytrocic index n increases. 

— For polynomial laws, the values of C/tiitR 2 increases from 0.24 to 0.4 when the power a increases 
and for the polytropic laws, the values of C/ttitR 2 decreases from 0.4 to 0 when the polytropic 
index n increases. 
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Online Supplement 
A Gravitational potentials 

A.l Potential of a homogeneous ellipsoid at an internal point 

Let us consider the contribution by the fc-th ellipsoid at one point = t^ix + r]ey + C,tz on the surface of the 
t -th ellipsoid, assumed interior to it. The potential of these ellipsoid may be written as 


(A.l) 


Vk(ft, fie, 0 ) — Kq k + A e £f + A v r/i + , 
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where the coefficients Vok, A^, A v and are 

2 r ° 

Vok — 7T G(7 k C k 

A.£ = 7T GfJkOtk / 

Jo 
r 

Arj — 7rG(T k i3 k 

A c = 7rGcrfc7fc 


dt 


o \Z0~+ a k t)( 1 + /3fcl)(l + f) 

dt 


(1 + CXkt)\J (1 + akt){ 1 + /3fct)(l+t) 

Z 100 _dt_ 

o (1 + /3fct) \/ (1 + cxkt) (1 + /3fct) (1 + t) 

r°° _ dt _ 

0 (1 + 7kt)\/0- + a k t)( 1 + 0kt) (1 + t) 


(A.2) 


and G is the gravitational constant (see Tisserand, 1891, Chap. 8 and 13; Jardetzky, 1958, Sec. 2.2). Then the 
derivatives of the potential are 


1 dV k _ 

r>.j- — 27T (jrO'kQ'k 

& oft 


dt 


1 dV k 
r/e dpt 


— 2'ivGo'k[3k 


= 2nGa k 

(e dCe 


l 0 (l + a fc t) 3 /2(i + /3 ^)i/2(i +t )i/2 

r°° dt 

1 o (l + a k ty/Z(l + p k t) 3/2 (l + t) 1/2 
0 dt 


I o (l + a fc f) 1 /2( 1 +/3 fc t)i/2(l+f)3/2’ 

and its contribution to the t- th equation of equilibrium are 
(i) nr \ 1 dVk at dVk 




— 27vGak 


dt. 


ctk(l - ae)t + (ctk - at) 

Jo (l + a fc t) 3 / 2 (l + /? fc t) 1 /2(i + i )3/2 
Neglecting terms of order 2 in the flattenings we can write 

X^(Vfc) = 2nGa k \(l - ae)Si(a k , /3k) + ( a k - at)S 2 (a k , /3k) 


where Si are functions of ak and f3k- If we expand around ak,(3k ~ 1 

r°° t 

Si(ak,Pk) = 


I o (1 wW+fttW+t) 3 / 2 

- 4 + l (1_afc)+ l55 (1_/3fc) + --- 


Si{a k , j3 k ) = 


'o (1 + a fc t) 3 / 2 (l + / 5 fc 0 1/2 (l + 0 3/2 


dt 


dt 




Finally, we obtain 


If k = t we have 


and similar equation for xf? ■ 


(!) nr \ Gm k 
Xt \Vk) = -^3- 


2j6£ — 


6efc 




Gme 4q 


(A.3) 


(A.4) 


(A.5) 


(A.6) 

(A.7) 

(A.8) 


A.2 Potential of a homogeneous ellipsoid at an external point 

The potential generated by the fc-th homogeneous ellipsoid at an external point = fix + pty + Qz on the 
surface of the Atli layer may be presented by Laplace series. Neglecting harmonics of degree higher than 2 we 
have 


Vk(ft, rn, <U) = ^ [hid + J vvvl + kcCt + 2 I €v bw + 2I^pifi + 2I a f e ft] , (A.9) 
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where r* = |xf|, I k is the moment of inertia of the fc-th ellipsoid, relative to the center of mass and hj are 
the components of its inertia tensor (see Beutler, 2005; Murray and Dermott, 1999). If the reference axes are 
oriented following the principal axes of inertia, then = 0 if i ^ j. We may define Ak = IBk = I-q-q and 
Ck = la- Because 2I k = A k + B k + Ck 


Vk(&,r]e,Q) = ~ 


Gm k G[C k — Ak) 


-r - 


r* 2 r * 5 

The principal moments of inertia of a homogeneous ellipsoid are 

1 


G(C k - B k ) 

2 r * 5 


(3??! - r* 2 ) . 


A k = \m k (b 2 k + cl) = \m k cl ( + 1 

5 5 

B k = \m k {a\ + cl) = \m k cl ( — + 1 


Oik 


Ck = \m k {al + bl) = \m k c k (— + 
5 5 \a k 


1 2 

7 m k c k 


o k (dk 
o: k + Pk 


Ck, 


(A.10) 


(A.ll) 


so its subtraction can be approximated 


C k - A k = (1 ak) f k C k » tkCk » 7 mkRlck 
Ok + Pk 5 

Ck — Bk = ^~^ k l ak Ck » MfcCfe « ImkRl^k. 

Ok+ Pk 5 


(A. 12) 


therefore, the potential can be written as 


Vk(&, ye, Q) = — 


Gmk GmkRlek 


(36 2 - r* 2 ) 


GTTl k R k fA k 

5r * 5 


(3 % 2 - r* 2 ) 


(A.13) 


and its derivatives are 

1 8Vk Gmk 2Gm k Rl 
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(A.14) 


Then we obtain 


W( > }_dVk _ a^dVk_ = }_9Vk _ }_dVk 2ee 8V k 
* 1 k d& Q d(e £e d£e Q d(e Q dQ 

6Gm k R 2 2Gm k 

~ 5r*® £fc + 


(A.15) 


where we have neglected the higher-order products t k te and eep k - Finally, by making the approximation r* sc Re, 
we obtain 


xf\v k ) = ^ 


Gmk 


2t£ 


_ 6 £fe ( Rk V 
Re) 


and, in the same way as before 


(A.16) 


( 2 ) / T , x Gm k 
Xe ’(V k ) = -jp- 


2pe 


_ 6 p k f Rk V 
Re) 


(A.17) 


A.3 Tidal potential 


If r = rx is the position of the mass M, the tidal potential at a point = £ex + yey + C,ez on the surface of 
the l -tli layer is (Lambeck 1980) 

C'Mr* 2 

Vu d = - 5 — P 2 (? ■ He) (A.18) 
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where r* = |x^| S 2 Re and P 2 is the Legendre polynomial of degree two. The differential acceleration of this 
point is 

GM 


-VV ud = — 


therefore their derivatives are 


- ■ 


1 dV ud 
ie die 
1 dVud 
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1 avtu 
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Finally, the contribution of the tide in the equilibrium equations of the Ath layer is 


xPiVut) = 


1 dVud _ ae dVud. 

it d( e 


£,t die 

3 GM 


2 GM 

H Is - e A 


(A.19) 


(A.20) 


(A.21) 


and, similarly, 


xf\Vu d ) = 


1 dVtid _ Pt dVu d 
it die 


Vt 9ve 
GM 

—W- 


(A.22) 


However, we discard terms containing te and fie, because when we calculate the flattenings of each layer, 
they appear multiplied by a factor of the same order as e or /x, therefore we obtain 


X?\Vud) = - 


3 GM 


and 


xf ) (Vtid) = 0 . 


(A.23) 
(A.24) 


B MacLaurin, Jeans and Roche ellipsoids 

In this appendix, we show that, in the homogeneous case, the flattenings defined by eqs. (1131) correspond to 
the MacLaurin and Jeans ellipsoids, respectively. In the homogeneous case, the equation ([9]) takes the simple 
form 

xiV) = X^(V tid )+x^(Vi) = 12 2 , (B.l) 

In this case, the contribution of the attraction of this ellipsoid to the equilibrium equation is 


and the tidal contribution is 


xF^i) 

X?Vi) 


Gm 4ei 
Gm 4/ki 


X ( e\Vu d ) 

xt\Vud) 


3GM 

~ r 3 ’ 
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The solutions can readily be found. They are 
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(B.4) 
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When the effect of the rotation is the only considered, the homogeneous solution was discovered by Newton 
and later generalized by MacLaurin (Chandrasekhar, 1969). The solution is an oblate spheroid with semiaxes 
a = b > c, with polar flattening ei = /ri = cm 


_ 5R 3 C 2 

4 mG 


(B.5) 


On the other hand, when the effect of the tide is the only considered, the homogeneous solution to the first 
order is a Roche ellipsoid with two equal axes (Tisserand, 1981). This equlibrium ellipsoid was later called 
Jeans ellipsoid (Chandrasekhar, 1969). It is a prolate spheroid with semiaxes a > b = c, with p,\ = 0 and 


ei = tj 


- 11M. (lY 

tJ 4 m \ r ) 


(B.6) 


Figure [9] illustrates the two cases. 

When we consider both effects, the solutions (IB.41) correspond to a triaxial ellipsoid with semiaxes a > b > c 
and the polar flattenings are ei = tj + £m and /ii = cm- 




Fig. [9ja Fig. [9jb 

Fig. 9 a) MacLaurin oblate spheroid when the axis of rotation is on the z axis. In this case, a = b > c. b) 
Jeans prolate spheroid when the mass disturbing the particle lies on the axis x. In this case, a > b = c. 


C Extension to the continuous case 


In this appendix, we show the intermediary cal cula tions necessary to find the Clairaut’s equation, following 
Tisserand (1891). Let us start from equations (11511 . and let us make the transition to the continuum. If we 
introduce the notation A(g k ) = gk — Pk- i and the boundary values Ro = 0 and p n +i = 0, we may rewrite the 
terms on the right hand side of eqn. m as 
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We also have 
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(B.3) 


In order to extend to the continuous case, we assume that the number of layers tends to infinity so that the 
increments Ag are infinitesimal quantities. When Ag —» 0, all terms in which the infinitesimals appear alone 
tends to zero. Only the terms in which they are integrated over an infinite number of layers remain. Hence, 


aekRk + ^ PtkRk 

k =1 k=£+1 

'll 


2tt 

m,T 

5m(R) 
2 m,T 


h L-o p{s)d ( s5w(s) ) + r3 f " 
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(B.4) 


where R„ is the mean radius of the ellipsoid and m(R) is the total mass enclosed by an ellipsoid of mean radius 
R. We thus have the integral equation 


7 (*>*(*> = i + S 


jp J s p{s)d(^s 5 'H{s)'j + R 3 ^ p(s)dH(s) 


(B.5) 


It is worth noting that the only restriction for the density profile is that it may be a piecewise continuous 
function. In order to normalize the problem, we first define a new variable 


R 

X ~ RJ 

so that *(0) = 0 in the center and x(R n ) = 1 on the surface. 

In the next step, we normalize the densities law 

p{R) = pop(x), 

where po is the density in the center and p(0) = 1. 

At last, we define the function 

f(x) = 3 [ p(z)z 2 dz, 

Jo 

where /(0) = 0 and /(1) = f n . The mass enclosed in the i^-sphere is then 

0 47T q 

m(R) = 4 n J p(r)r dr = —Rlpof{x), 


and the total mass is 


47T 3 

m-T = —Rnpojn- 


Then, the integral equation (lB~5l) becomes 

r 2 o f rz=x rz =1 

-f{x)'H{x) = ^-x 3 + p(z)d{z 5 H{z))+x 5 p{z)dH{z) 

^ J 2 = 0 J Z=X 

or, deriving with respect to * and dividing by 5x 4 , 

^R(x) + ^ + £ p(z)d'H(z). 

Deriving once more we obtain the Clairaut’s equation for the flattening profile R(x) 

6p(x)x 2 
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